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OPTIMIZED COMPUTATION WITH 
RECURSIVE POWER SERIES INTEGRATION 

ABSTRACT 

The power series technique is applied to the restricted 
three-body problem and recursion formulas are developed 
for the coefficients of the power series which reduce sig- 
nificantly the computation time required by the recursion 
formulas developed by previous authors. Van Flandern's 
method for maximizing the step size per unit of computer 
time is discussed and extended on the basis of relative error 
control. A method for controlling the magnitudes of the 
series coefficients is also presented. 



OPTIMIZED COMPUTATION WITH 
RECURSIVE POWER SERIES INTEGRATION 


INTRODUCTION 

The method of solution of second degree systems of differential equations 
by recurrent power series is well known [6], This technique is especially use- 
ful when variable step sizes are needed in the integration, and has the bonus 
property of automatically choosing the correct step size relative to some pre- 
set error tolerance. As pointed out by Fehlberg [2], a substantial saving in 
computing time is realizable under the above conditions by choosing a power 
series method over a Runge-Kutta or interpolation procedure. 

In this paper the power series technique will be applied to the restricted 
three-body problem and a few "tricks" will be presented which reduce signifi- 
cantly the computation time required by the recursion formulas presented in the 
papers we have seen [2, 5, 6]. Van Flandern’s method [7] for maximizing the 
step size per unit of computer time is discussed and extended on the basis of 
relative error control. A method for controlling the magnitudes of the series 
coefficients is also presented. 


THE RESTRICTED THREE-BODY PROBLEM 


Second degree systems of differential equations are not always met in prac- 
tice. However, upon introducing appropriate auxiliary functions a given system 
may in many instances be transformed into one of the second degree. It is to 
be noted that in general for a particular set of equations there are no unique 
auxiliary functions which will transform the system to the above form, but rather 
one must search for the functional substitutions which lead to more symmetric 
recursion formulas and fewer computational operations. As an illustration, con- 
sider the restricted problem of three bodies. 


In a rotating coordinate system the equations of motion are 
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where \i is the mass of the lighter body and 1 -j± is the mass of the heavier 
body. The commonly used auxiliary functions 
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upon introduction into the equations of motion yield the following second degree 
system: 
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and substituting into (3) we obtain, after equating coefficients of like powers, 
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The initial conditions of the problem yield x Q , x 1S y 0 , and y x , and r Q , s 0 , u 0 , 
v 0 are obtained from (2), the auxiliary functions. It may be easily seen that 
proceeding through this algorithan to calculate the (n + 1) st term, i.e., x n + t 
and y n + j , there are required 

14n + 4 additions, 

16n + 16 multiplications, 
and 6 divisions. 


Further, 26n + 4 quantities are addressed. 

Let us return now to (1) and consider instead the following auxiliary func- 
tions: 
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Substituting (6) into (1) then yields the second degree system: 
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and substituting into (7) there results 
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for n>l. Again, the initial conditions yield s 0 , s x , y 0 , and y x ; and a Q , fi Q , p Q , 
q 0 are determined from (6). However, the calculation of the (n + 1) st term now 
involves 


9n + 5 additions 

7n + 20 multiplications 

and 4 divisions. 

Also, only 12n + 12 quantities are addressed. Thus, by simply choosing the 
auxiliary functions (6) instead of (2) the computation time for high order terms 
will be approximately halved. It is to be noted that symmetries such as in the 
first of equations (9) occur frequently in practice and full advantage of these 
symmetries should be taken in a computational scheme. 
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MINIMIZATION OF COMPUTER TIME 


In a power series integration method, the solution is represented as a power 
series, e.g., „ 



n = 0 


For a given number of terms n, t is chosen so that x has a desired accuracy. By 
using the solutions of the previous step as initial conditions for a new step and 
computing a new power series, the solution is extended by analytic continuation. 
The following technique, due to Thomas Van Flandern,[ 7 ] minimizes the computer 
time required to integrate over a given interval. Let r designate total compu- 
tation time for a step and assume that a fixed time, r 0 , is the initialization time 
for each step. Since the computation time required for the (n + 1) st. term is 
greater than that for the n th term, r - t (n) is an increasing monotone function. 
Further, with the addition of more terms in the series the integration step size 
h increases for fixed accuracy and approaches h max , the radius of convergence. 
Thus h = h(n) is also an increasing monotone function. 

We wish to make h as large as possible for a given amount of computer time. 
Thus we are led to choose the value of n such that the quantity h/r is a maxi- 
mum. That this ratio has at least one maximum is obvious from the fact that 
h/r is zero for n = 0, approaches zero as ngets large, and has only positive values 
for n> 0. 

The ratio h/r should be calculated for each new term of the series, and a 
check made for its maximum value. The step size h may be estimated from the 
newly calculated term and a preset error tolerance, while r may be obtained 
from an accurate computer clock or from operation counting for the particular 
recursive formulas used. 

The following graphs illustrate the above procedure: 
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In practice, convergence of the series will usually insure that the sum of 
neglected terms will be smaller in magnitude than the first term omitted. In 
such cases we may calculate the step size for each new term from 


h = 



( 11 ) 


where a is the n th coefficient of the series and e is the absolute error tolerance. 

n 

Many times with this procedure it is found to be true that the magnitudes of the 
calculated series coefficients are such that h = h(n) is not an increasing mono- 
tone function, i.e., it has small oscillations in violation of the criteria of appli- 
cability for Van Flandern's procedure. However, the authors have found in gen- 
eral that such coefficients which disrupt the smoothly increasing sequence of the 
magnitude of h are isolated and so the requirement that the ratio h/r of three 
consecutive terms form a strictly decreasing sequence yields excellent results 
as a test for the maximum value of h/r. 


In order to avoid unnecessary labor in computation it is desirable that all 
errors committed at any integration step contribute uniformly to the entire in- 
tegration error. Thus the error tolerance e of equation (11) must be a relative, 
rather than an absolute error. The authors have found it convenient to form the 
absolute error, e a , for each integration step by selecting a relative error e f for 
the entire integration and from the initial conditions a 0 at each step, forming 
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When a 0 is small, however, e a is replaced by a preset value. 

A major problem which occurs in practice is that the series coefficients a 
may become quite large for increasing n. This results in a loss of accuracy when 
the truncated series is evaluated at a value of the independent variable for a given 
integration step. However, the following convenient property of recursive power 
series may be used to easily overcome this difficulty. 


Given the series (10), i.e. 
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we may scale t such that t = hs so that the above series becomes 


CO 
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where the coefficients b are calculated with the recursive relations which 

n 

determine a slightly changed and with simple changes in the initial conditions 
f 1 ! . M. By assuming that the step size for a given integration step will be 
nearly equal to that of the previous step, we set s equal to the previous step 
size so that h ~ 1 and the coefficients b n are well behaved. 

In the interest of reducing roundoff error over the course of the integration 
the authors have found it useful to set the scale factor equal to a power of two 
so that the operations involving the scale factor are merely shifts within a bi- 
nary digital computer. 

We have programmed the above techniques for the restricted problem of 
three bodies and applied the same initial conditions as those of Fehlberg [ 3 ]for 
the integration of a periodic orbit. Our inability to obtain a replica of the opera- 
ting program used by this author has led to limited comparison. However, 
our own results indicate that the reduced computational procedures outlined in 
this paper produce a 40%-50% time saving over a conventional power series 
integration of this periodic orbit with each step constrained to sixteen terms, 
where both methods were performed in double precision arithmetic on an IBM 
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7094 computer with a preset error tolerance of 10~ 17 . We conclude from this that 
the above techniques lead to the same high precision in comparable computing 
time as the high order Runge-Kutta tranformation type formulas of Fehlbergt 3 ! 
with the advantage of simpler and more compact programming. 
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